/*********** NOTES HEADER *********** 

DESCRIPTION: This program estimates the logit
	and the marginal effects separately by merger

Datasets used:
 (1) MatchedSample_Full

Datasets created:
 (1) MbM_hyp_zip. 
 (2) MbM_diab_zip. 
 (3) MbM_death_zip 
 
************************************/ 

cd "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Program\PropensityScore\Patients\ResponsetoReferee\OtherPropScoreMethods\BiasAdjustment"
set more off
capture log close
set matsize 10000
clear all

global origData "N:\MedicareClaims-P045601-BE"
global dataIn "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Data-In\"
global dataOut "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Data-Out\"
global dpath "N:\MedicareClaims-P045601-BE\Work\ay_data"
global dataProp "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Data-Out\PropScore\Patients"
global logs "N:\MedicareClaims-P045601-BE\Work\hosp_retro\health_out\Logs\PropScore\Patients"
global dpath "N:\MedicareClaims-P045601-BE\Work\ay_data"
global skapath "N:\MedicareClaims-P045601-BE\Work\ska"

adopath +  "N:/SIL-Common/estout"
adopath +  "N:/SIL-Common/outreg2"
adopath +  "N:/SIL-Common/reghdfe-master/package"

log using "./5.PS_mergerbymerger.log", replace

program getColnames 
 quietly {
  matrix A = e(b)'
 
 local rownames : rowfullnames A
 local c : word count `rownames'
 local names : colfullnames A

 svmat A, names(col)

 gen rownames = ""
 forvalues i = 1/`c' {

  replace rownames = "`:word `i' of `rownames''" in `i' 
   } 

 }

end


program getRegStats
 quietly {
  matrix A = e(b)'
 
 local rownames : rowfullnames A
 local c : word count `rownames'
 local n = `c' + 1
 local pll = `c' + 2
 local rsq = `c' + 3
 local ndrop = `c' + 4

 replace y1 =e(N) in `n' 
 replace y1 =e(ll) in `pll' 
 replace y1 =e(r2_p) in `rsq' 
 replace y1 =e(N_drop) in `ndrop' 

 }

end

use ./MatchedSample_Full

	*Use aggregate
	capture drop VMP0
	gen VMP0 =0
	levelsof vm_provider1 if inlist(vm_provider1, 1,12,20,22,34,37 ), local(vmp)
	foreach i of local vmp {
	replace VMP0 =1 if vm_provider1==`i'
	replace VMP0 =1 if vm_provider2==`i'
	replace VMP0 =1 if vm_provider3==`i'
	}

	capture drop pmVMP0
	gen pmVMP0 = VMP0*postmerger

	replace vm_provider1=0 if vm_provider1==.

	levelsof vm_provider1 if !inlist(vm_provider1, 1,12,20,22,34,37 ), local(vmp)



*Combined sexes
		foreach var of varlist oc_glaucoma oc_diabcomp1 oc_diabcomp2 {
		*This limits the sample to diabetes patients
		 forval j=3(1)3  {

		if regexm("`var'","oc_([a-z0-9]+)$") { 
		  local pre = "`=regexs(1)'"

		 capture drop zip3`pre' zipct`pre'
		 gen zip3`pre'=zip3dig 
		 egen zipct`pre'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
		 replace zip3`pre'=0 if zipct`pre'<1250
		 replace zip3`pre'=0 if kv_`pre'==0 
		 replace zip3`pre'=STATE_encode+1000 if zip3`pre'==0 
		 
		 capture rename mxicd3_chronic mxicd3_endocrine 

		logit `var' pm* VMP* vmerger##c.nyear i.vmerger##i.sex_pat i.vmerger#c.nyear nonwhite##i.agecat nonwhite#i.sex_pat i.p_age i.agecat##i.sex_pat i.qtr  ib452.zip3`pre' fmp pspcl30 pspcl06 pspcl69 pspcl18 pspcl48 pspcl93 pspclOther mxicd1_acute mxicd5_acute mxicd18_acute mxicd*_c* fs* lnvisit* [pweight=samp_weight] if flag`j'==1 & kv_`pre'==1 , cluster( propScoreId )noomit

			foreach i of local vmp {
				nlcom exp(_cons+_b[pmVMP`i']+_b[VMP`i'] + _b[77.p_age] )/(1+exp(_cons+_b[pmVMP`i']+_b[VMP`i'] + _b[77.p_age] )) - (exp( _cons+ _b[77.p_age] +_b[VMP`i'])/(1+exp( _cons+_b[77.p_age] + _b[VMP`i'] )))
				  mat A`i'=r(b)
				  mat B`i'=r(V)
				  mat C`i'=[A`i',B`i']
				  
				   if `i'==0 {
				   mat D=C`i''
				   }
				   
				   if `i'!=0 {
				   mat D=[D,C`i'']
				   }
			  }


		  matrix beta  = e(b)' 
		  matrix sesq = vecdiag(e(V))' 
		  svmat beta
		  getColnames
		  getRegStats
		  svmat sesq
		  rename rownames rn`pre'_`j'
		  capture drop y1 
		  capture drop j 
		  capture drop i
		  rename beta1 beta_`pre'_`j'
		  rename sesq1 sesq_`pre'_`j'
		  matrix Dt = D'
		  svmat Dt
		  rename Dt1 mem_beta_`pre'_`j'
		  rename Dt2 mem_sesq_`pre'_`j'
		   }
		  }
		}

		save "./MbM_diab_zip.dta", replace


	rename mxicd3_endocrine mxicd3_chronic 

		foreach var of varlist oc_acuteicd8 oc_ihd oc_ami {
		*This limits the sample to hyptertensive patients
		 forval j=2(1)2  {

		if regexm("`var'","oc_([a-z0-9]+)$") { 
		  local pre = "`=regexs(1)'"

		 capture drop zip3`pre' zipct`pre'
		 gen zip3`pre'=zip3dig 
		 egen zipct`pre'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
		 replace zip3`pre'=0 if zipct`pre'<1250
		 replace zip3`pre'=0 if kv_`pre'==0 
		 replace zip3`pre'=STATE_encode+1000 if zip3`pre'==0 
		 
		 capture rename mxicd8_chronic mxicd3_circulatory 

		*levelsof idpost if idpost!=0, local(idlocal) 

		logit `var' pm* VMP* vmerger##c.nyear i.vmerger##i.sex_pat i.vmerger#c.nyear nonwhite##i.agecat nonwhite#i.sex_pat i.p_age i.agecat##i.sex_pat i.qtr  ib452.zip3`pre' fmp pspcl30 pspcl06 pspcl69 pspcl18 pspcl48 pspcl93 pspclOther mxicd1_acute mxicd5_acute mxicd18_acute mxicd*_c* fs* lnvisit* [pweight=samp_weight] if flag`j'==1 & kv_`pre'==1 , cluster( propScoreId ) noomit

			foreach i of local vmp {
				nlcom exp(_cons+_b[pmVMP`i']+_b[VMP`i'] + _b[77.p_age] )/(1+exp(_cons+_b[pmVMP`i']+_b[VMP`i'] + _b[77.p_age] )) - (exp( _cons+ _b[77.p_age] +_b[VMP`i'])/(1+exp( _cons+_b[77.p_age] + _b[VMP`i'] )))
				  mat A`i'=r(b)
				  mat B`i'=r(V)
				  mat C`i'=[A`i',B`i']
				  
				   if `i'==0 {
				   mat D=C`i''
				   }
				   
				   if `i'!=0 {
				   mat D=[D,C`i'']
				   }
			  }


		  matrix beta  = e(b)' 
		  matrix sesq = vecdiag(e(V))' 
		  svmat beta
		  getColnames
		  getRegStats
		  svmat sesq
		  rename rownames rn`pre'_`i'
		  capture drop y1 
		  capture drop j 
		  capture drop i
		  rename beta1 beta_`pre'_`i'
		  rename sesq1 sesq_`pre'_`i'
		  matrix Dt = D'
		  svmat Dt
		  rename Dt1 mem_beta_`pre'_`i'
		  rename Dt2 mem_sesq_`pre'_`i'
		 
		   }
		  }
		}

		save "./MbM_hyp_zip.dta", replace

	rename  mxicd8_circulatory  mxicd8_chronic 

 
		foreach var of varlist oc_death   {
		 forval j=1(1)3  {


		if regexm("`var'","oc_([a-z0-9]+)$") { 
		  local pre = "`=regexs(1)'"

		 capture drop zip3`pre' zipct`pre'
		 gen zip3`pre'=zip3dig 
		 egen zipct`pre'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
		 replace zip3`pre'=0 if zipct`pre'<1500
		 replace zip3`pre'=0 if kv_`pre'==0 
		 replace zip3`pre'=STATE_encode+1000 if zip3`pre'==0 
		 
		 drop zipct`pre'
		 egen zipct`pre'=count(zip3dig) if kv_`pre'==1 , by(zip3dig)
		 replace zip3`pre'=2000 if zipct`pre'<1250
		  
		logit `var' pm* VMP* vmerger##c.nyear i.vmerger##i.sex_pat i.vmerger#c.nyear nonwhite##i.agecat nonwhite#i.sex_pat i.p_age i.agecat##i.sex_pat i.qtr  ib452.zip3`pre' fmp pspcl30 pspcl06 pspcl69 pspcl18 pspcl48 pspcl93 pspclOther mxicd1_acute mxicd5_acute mxicd18_acute mxicd*_c* fs* lnvisit* [pweight=samp_weight] if flag`j'==1 & kv_`pre'==1 , cluster( propScoreId ) noomit

			foreach i of local vmp {
				nlcom exp(_cons+_b[pmVMP`i']+_b[VMP`i'] + _b[77.p_age] )/(1+exp(_cons+_b[pmVMP`i']+_b[VMP`i'] + _b[77.p_age] )) - (exp( _cons+ _b[77.p_age] +_b[VMP`i'])/(1+exp( _cons+_b[77.p_age] + _b[VMP`i'] )))
				  mat A`i'=r(b)
				  mat B`i'=r(V)
				  mat C`i'=[A`i',B`i']
				  
				   if `i'==0 {
				   mat D=C`i''
				   }
				   
				   if `i'!=0 {
				   mat D=[D,C`i'']
				   }
			  }



		  matrix beta  = e(b)' 
		  matrix sesq = vecdiag(e(V))' 
		  svmat beta
		  getColnames
		  getRegStats
		  svmat sesq
		  rename rownames rn`pre'_`j'
		  capture drop y1 
		  capture drop j 
		  capture drop i
		  rename beta1 beta_`pre'_`j'
		  rename sesq1 sesq_`pre'_`j'
		  matrix Dt = D'
		  svmat Dt
		  rename Dt1 mem_beta_`pre'_`j'
		  rename Dt2 mem_sesq_`pre'_`j'
		 
		   }
		  }
		}

		save "./MbM_death_zip.dta", replace
		clear
 



